#include "solver_tpw.h"

// QtCore
#include <QtCore/QDir>
#include <QtCore/QFile>
#include <QtCore/QTextCodec>
#include <QtCore/QDataStream>
#include <QtCore/QStringBuilder>
#include <QtCore/QCoreApplication>
// Локальные
#include "settings.h"
#include "math/math.h"
#include "linal/los.h"
#include "linal/bcg.h"
#include "solver_tpwt.h"
#include "fem/trianglefem.h"
#include "fem/triangulation.h"
#include "geometry/2d/vector.h"

//================================================================================================================================
//================================================================================================================================
SolverTPW::SolverTPW() : m_NetSearcher(NULL)
{
    //****************************************************************************************************************************
    // Создание и очистка директории для вывода результатов.
    //****************************************************************************************************************************
    m_OutputPath = Settings::instance()->setting("OutputPath").toString();
    if(!QDir::isAbsolutePath(m_OutputPath))
        m_OutputPath = QDir::cleanPath(m_OutputPath).prepend(qApp->applicationDirPath() + "/");

    if(!QDir(m_OutputPath).exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputPath);
    // TODO: Очистка директории.

    //****************************************************************************************************************************
    // Считывание физических и геометрических параметров задачи. Перевод скорости из см/ч в м/с.
    //****************************************************************************************************************************
    m_TMax = Settings::instance()->setting("Physics/MaxTemperature").toDouble();
    m_TMin = Settings::instance()->setting("Physics/MinTemperature").toDouble();

    //****************************************************************************************************************************
    // Считывание параметров решателей.
    //****************************************************************************************************************************
    m_MaxiterSLAE = Settings::instance()->setting("Math/SlaeIters").toInt();
    m_MaxiterNonlin = Settings::instance()->setting("Math/NonlinIters").toInt();

    m_EpsSLAE = Settings::instance()->setting("Math/SlaeResidual").toDouble();
    m_EpsNonlin = Settings::instance()->setting("Math/NonlinResidual").toDouble();

    //****************************************************************************************************************************
    // Инициализация решателей.
    //****************************************************************************************************************************
    m_Solvers.append(new linal::LOS(m_MaxiterSLAE, m_EpsSLAE));
    m_Solvers.append(new linal::BCG(m_MaxiterSLAE, m_EpsSLAE));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_DiffStep = Settings::instance()->setting("Math/DiffStep").toDouble();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_W = Settings::instance()->setting("Math/RelaxFactor").toDouble();

    //****************************************************************************************************************************
    // Инициализация поисковика по сетке.
    //****************************************************************************************************************************
    m_NetSearcher = new fem::TriangleNetSearch(60, 60, 0.0025, 0.0025);

    //****************************************************************************************************************************
    // Создание конечноэлементной сетки.
    //****************************************************************************************************************************
    fem::Triangulation triangulation(QStringList() << "T" << "Psi" << "Omega");

    int basis_points[4];
    basis_points[0] =  triangulation.addBasisPoint(geometry::_2d::Point(0.05, 0.00));
    basis_points[1] =  triangulation.addBasisPoint(geometry::_2d::Point(0.10, 0.00));
    basis_points[2] =  triangulation.addBasisPoint(geometry::_2d::Point(0.05, 0.05));
    basis_points[3] =  triangulation.addBasisPoint(geometry::_2d::Point(0.10, 0.05));

    int border_segments[4];
    border_segments[0] =  triangulation.addBorderSegment(basis_points[0], basis_points[1], 99);
    border_segments[1] =  triangulation.addBorderSegment(basis_points[0], basis_points[2], 99);
    border_segments[2] =  triangulation.addBorderSegment(basis_points[1], basis_points[3], 99);
    border_segments[3] =  triangulation.addBorderSegment(basis_points[2], basis_points[3], 99);

    triangulation.addBorderSegmentBCInfo("T", border_segments[1], 1, 0);
    triangulation.addBorderSegmentBCInfo("T", border_segments[2], 1, 1);

    triangulation.addBorderSegmentBCInfo("Psi", border_segments[0], 1, 0);
    triangulation.addBorderSegmentBCInfo("Psi", border_segments[1], 1, 0);
    triangulation.addBorderSegmentBCInfo("Psi", border_segments[2], 1, 0);
    triangulation.addBorderSegmentBCInfo("Psi", border_segments[3], 1, 0);

    triangulation.addBorderSegmentBCInfo("Omega", border_segments[0], 1, 0);
    triangulation.addBorderSegmentBCInfo("Omega", border_segments[1], 1, 1);
    triangulation.addBorderSegmentBCInfo("Omega", border_segments[2], 1, 2);
    triangulation.addBorderSegmentBCInfo("Omega", border_segments[3], 1, 3);

    triangulation.addQuadrangleSubdomain(border_segments[3],  border_segments[2],  border_segments[0],  border_segments[1], 0);

    m_Net = triangulation.triangulate();
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //Формирование информации о краевых условиях.
    //****************************************************************************************************************************
    triangulation.getDirichletBoundaryInfo("T", m_TopTBC, m_TypeTBC);
    triangulation.getDirichletBoundaryInfo("Psi", m_TopPsiBC, m_TypePsiBC);
    triangulation.getDirichletBoundaryInfo("Omega", m_TopOmegaBC, m_TypeOmegaBC);

    //****************************************************************************************************************************
    // Формирование начальных значений для векторов решения.
    //****************************************************************************************************************************
    m_T = linal::Vector(m_Net.numOfPoints(), m_TMax);
    m_TOldNolin = linal::Vector(m_Net.numOfPoints(), m_TMax);

    m_Vr = linal::Vector(m_Net.numOfPoints(), 0);
    m_Vz = linal::Vector(m_Net.numOfPoints(), 0);

    m_Psi = linal::Vector(m_Net.numOfPoints(), 0);
    m_PsiOldNolin = linal::Vector(m_Net.numOfPoints(), 0);

    m_Omega = linal::Vector(m_Net.numOfPoints(), 0);
    m_OmegaOldNolin = linal::Vector(m_Net.numOfPoints(), 0);

    //****************************************************************************************************************************
    // Формирование портрета СЛАУ. Выделение памяти под СЛАУ.
    //****************************************************************************************************************************
    m_F = linal::Vector(m_Net.numOfPoints());

    QVector<int> igl, jgl, igu, nju;
    fem::triangles::formMatrixProfile(m_Net, igl, jgl, igu, nju);
    m_A = linal::CSRMatrix(m_Net.numOfPoints(), igl, jgl, igu, nju);

    foreach(linal::SLAESolver* solver, m_Solvers)
        solver->adaptForSLAE(m_A.dimension());

    //****************************************************************************************************************************
    // Расплав.
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Liquid));
    m_Materials[0].setProperty(Material::Diffusivity, 0);
    m_Materials[0].setProperty(Material::Density, Settings::instance()->setting("Liquid/Density").toDouble());
    m_Materials[0].setProperty(Material::MolarMass, Settings::instance()->setting("Liquid/MolarMass").toDouble());
    m_Materials[0].setProperty(Material::Viscosity, Settings::instance()->setting("Liquid/Viscosity").toDouble());
    m_Materials[0].setProperty(Material::HeatCapacity, Settings::instance()->setting("Liquid/HeatCapacity").toDouble());
    m_Materials[0].setProperty(Material::Expansiveness, Settings::instance()->setting("Liquid/Expansiveness").toDouble());
    m_Materials[0].setProperty(Material::HeatConductivity, Settings::instance()->setting("Liquid/HeatConductivity").toDouble());
    m_Materials[0].setProperty(Material::CrystallizationHeat, Settings::instance()->setting("Liquid/CrystallizationHeat").toDouble());
}

//================================================================================================================================
//================================================================================================================================
bool SolverTPW::solveSLAE(const linal::CSRMatrix &A, const linal::Vector &f, linal::Vector &x)
{
    linal::Vector x_old = x;

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveLU(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveDi(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solve(A,f,x) < m_EpsSLAE)
            return true;
    }

    qWarning("SLAE wasn't solved!");
    return false;
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::solve()
{
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int nolin_iter(0);
    double nolin_eps(0.);

    do
    {
        formV();

        //Решаем уравнение для T
        m_A.zeroize();
        m_F.zeroize();
        formTSLAE();
        solveSLAE(m_A, m_F, m_T);
        m_T = (1 - m_W)*m_TOldNolin + m_W*m_T;

        //Решаем уравнение для Omega
        m_A.zeroize();
        m_F.zeroize();
        formOmegaSLAE();
        solveSLAE(m_A, m_F, m_Omega);
        m_Omega = (1 - m_W)*m_OmegaOldNolin + m_W*m_Omega;

        //Решаем уравнение для Psi
        m_A.zeroize();
        m_F.zeroize();
        formPsiSLAE();
        solveSLAE(m_A, m_F, m_Psi);
        m_Psi = (1 - m_W)*m_PsiOldNolin + m_W*m_Psi;

        //Рассчитываем невязку по нелинейности
        double res_t = linal::Vector::normSquared(m_T - m_TOldNolin)/linal::Vector::normSquared(m_T);
        double res_psi = linal::Vector::normSquared(m_Psi - m_PsiOldNolin)/linal::Vector::normSquared(m_Psi);
        double res_omega = linal::Vector::normSquared(m_Omega - m_OmegaOldNolin)/linal::Vector::normSquared(m_Omega);
        nolin_eps = qMax(qSqrt(res_t)/m_W, qMax(qSqrt(res_psi)/m_W, qSqrt(res_omega)/m_W));

        m_TOldNolin = m_T;
        m_PsiOldNolin = m_Psi;
        m_OmegaOldNolin = m_Omega;

        cout << nolin_iter << ": " << nolin_eps << endl;
    }
    while(++nolin_iter < m_MaxiterNonlin && nolin_eps > m_EpsNonlin);
    writeResult();
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::formV()
{
    m_Vr.resize(m_Net.numOfPoints());
    m_Vz.resize(m_Net.numOfPoints());

    //Находим r-компоненту скорости
    for(int i = 0; i < m_Vr.dimension(); i++)
    {
        if(!m_TopPsiBC.contains(i))
        {
            geometry::_2d::Point p(m_Net.point(i).shifted(0, -m_DiffStep));
            int triangle = m_NetSearcher->findTriangle(m_Net, p);
            m_Vr[i] = - (fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/(-m_DiffStep);
        }
        else
            m_Vr[i] = 0;
    }

    //Находим z-компоненту скорости
    for(int i = 0; i < m_Vz.dimension(); i++)
    {
        if(!m_TopPsiBC.contains(i))
        {
            geometry::_2d::Point p(m_Net.point(i).shifted(m_DiffStep, 0));
            int triangle = m_NetSearcher->findTriangle(m_Net, p);

            if(m_Net.point(i).x() == 0)
                m_Vz[i] = 2.*(fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/m_DiffStep;
            else
            {
                m_Vz[i] = (fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/m_DiffStep;
                m_Vz[i] += m_Psi[i]/m_Net.point(i).x();
            }
        }
        else
            m_Vz[i] = 0;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::tBC()
{
    for(int i = 0; i < m_TopTBC.size(); i++)
    {
        m_A.makeIdentity(m_TopTBC[i]);

        switch(m_TypeTBC[i])
        {
        case 0:
            m_F[m_TopTBC[i]] = m_TMax;
            break;
        case 1:
            m_F[m_TopTBC[i]] = m_TMin;
            break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::psiBC()
{
    for(int i = 0; i < m_TopPsiBC.size(); i++)
    {
        m_A.makeIdentity(m_TopPsiBC[i]);
        m_F[m_TopPsiBC[i]] = 0;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::omegaBC()
{
    for(int i = 0; i < m_TopOmegaBC.size(); i++)
    {
        m_A.makeIdentity(m_TopOmegaBC[i]);

        switch(m_TypeOmegaBC[i])
        {
        case 0: // Дно
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x();
            double z = m_Net.point(m_TopOmegaBC[i]).y() + m_DiffStep;

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        case 1: // Левая симметрии
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x() + m_DiffStep;
            double z = m_Net.point(m_TopOmegaBC[i]).y();

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        case 2: // Правая стенка
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x() - m_DiffStep;
            double z = m_Net.point(m_TopOmegaBC[i]).y();

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        case 3: // Верх
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x();
            double z = m_Net.point(m_TopOmegaBC[i]).y() - m_DiffStep;

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
            break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::formTSLAE()
{
    double density;
    double heat_capacity;
    double heat_conductivity;
    double D[3][3];             // Матрица, используемая для определения базисных функций
    double C[3][3][3];          // Локальная матрица массы
    double detD;                // Определитель матрицы D
    double abs_detD;            // Модуль определетиля матрицы D
    double r[3];                // r-координаты вершин КЭ
    double z[3];                // z-координаты вершин КЭ
    int top[3];                 // Индексы вершин КЭ

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        density = m_Materials[m_Net.material(cur_triangle)].property(Material::Density);
        heat_capacity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatCapacity);
        heat_conductivity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatConductivity);

        for(int i = 0; i < 3; i++)
        {
            r[i] = m_Net.point(cur_triangle, i).x();
            z[i] = m_Net.point(cur_triangle, i).y();
            top[i] = m_Net.pointIndex(cur_triangle, i);
        }

        detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
        abs_detD = qAbs(detD);

        D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
        D[0][1] = (z[1] - z[2])/detD;
        D[0][2] = (r[2] - r[1])/detD;
        D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
        D[1][1] = (z[2] - z[0])/detD;
        D[1][2] = (r[0] - r[2])/detD;
        D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
        D[2][1] = (z[0] - z[1])/detD;
        D[2][2] = (r[1] - r[0])/detD;

        C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
        C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
        C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

        C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
        C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
        C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

        C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
        C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
        C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                double addition(0);
                addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*heat_conductivity*abs_detD/6.;
                addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                             m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                             m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1]*heat_capacity*density;
                addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                             m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                             m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2]*heat_capacity*density;
                m_A(top[i], top[j]) += addition;
            }
        }
    }

    //****************************************************************************************************************************
    //Учитываем краевые условия.
    //****************************************************************************************************************************
    tBC();
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::formPsiSLAE()
{
    double D[3][3];         // Матрица, используемая для определения базисных функций
    double C[3][3][3];      // Локальная матрица массы
    double detD;            // Определитель матрицы D
    double abs_detD;        // Модуль определетиля матрицы D
    double r[3];            // r-координаты вершин КЭ
    double z[3];            // z-координаты вершин КЭ
    int top[3];             // Индексы вершин КЭ

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_A(top[i], top[j]) += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*abs_detD/6.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_F[top[i]] += m_Omega[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2]);

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_PsiOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])/r[j];
                    else
                    {
                        geometry::_2d::Point p(m_Net.point(top[j]).shifted(m_DiffStep, 0));
                        int triangle = m_NetSearcher->findTriangle(m_Net, p);
                        double psi_diff = (fem::triangles::functionValue(triangle, m_Net, m_PsiOldNolin, p) - m_PsiOldNolin[top[j]])/m_DiffStep;

                        m_F[top[i]] -= psi_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2]);
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    psiBC();
}

//================================================================================================================================
//================================================================================================================================
void SolverTPW::formOmegaSLAE()
{
    double viscosity;
    double expansiveness;
    double D[3][3];         // Матрица, используемая для определения базисных функций
    double C[3][3][3];      // Локальная матрица массы
    double detD;            // Определитель матрицы D
    double abs_detD;        // Модуль определетиля матрицы D
    double r[3];            // r-координаты вершин КЭ
    double z[3];            // z-координаты вершин КЭ
    int top[3];             // Индексы вершин КЭ

    //*******************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //*******************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            viscosity = m_Materials[m_Net.material(cur_triangle)].property(Material::Viscosity);
            expansiveness = m_Materials[m_Net.material(cur_triangle)].property(Material::Expansiveness);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*viscosity*abs_detD/6.;
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    addition -= m_Vr[top[0]]*C[i][j][0];
                    addition -= m_Vr[top[1]]*C[i][j][1];
                    addition -= m_Vr[top[2]]*C[i][j][2];

                    m_A(top[i], top[j]) += addition;
                }
            }

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] -= m_T[top[j]]*r[k]*(C[i][k][0] + C[i][k][1] + C[i][k][2])*D[j][1]*expansiveness*9.8;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_OmegaOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity/r[j];
                    else
                    {
                        geometry::_2d::Point p(m_Net.point(top[j]).shifted(m_DiffStep, 0));
                        int triangle = m_NetSearcher->findTriangle(m_Net, p);
                        double omega_diff = (fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, p) - m_OmegaOldNolin[top[j]])/(m_DiffStep);
                        m_F[top[i]] -= omega_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity;
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    omegaBC();
}

//================================================================================================================================
//Вывод результата в файл.
//================================================================================================================================
void SolverTPW::writeResult()
{
    QFile file;
    QTextStream ostream(&file);
    ostream.setRealNumberPrecision(12);

    //****************************************************************************************************************************
    //Выводим конечноэлементую сетку.
    //****************************************************************************************************************************
    file.setFileName(QString("/nvtr.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfElements(); ++i)
    {
        for(int j = 0; j < 3; ++j)
            ostream << m_Net.pointIndex(i, j) << " ";
        ostream << m_Net.material(i) << endl;
    }
    file.close();

    file.setFileName(QString("/rz.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); ++i)
        ostream << m_Net.point(i) << endl;
    file.close();

    //****************************************************************************************************************************
    //Выводим значения температуры.
    //****************************************************************************************************************************
    file.setFileName(QString("/T.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i) << " " << m_T[i] << endl;
    file.close();

    //************************************************************************************************************************
    //Выводим значения функции тока.
    //************************************************************************************************************************
    file.setFileName(QString("/Psi.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i) << " " << m_Psi[i] << endl;
    file.close();

    //************************************************************************************************************************
    //Выводим значения вихря.
    //************************************************************************************************************************
    file.setFileName(QString("/Omega.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i) << " " << m_Omega[i] << endl;
    file.close();

    //************************************************************************************************************************
    //Выводим значения r-компоненты скорости.
    //************************************************************************************************************************
    file.setFileName(QString("/Vr.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i) << " " << m_Vr[i] << endl;
    file.close();

    //************************************************************************************************************************
    //выводим значения z-компоненты скорости.
    //************************************************************************************************************************
    file.setFileName(QString("/Vz.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream << m_Net.point(i) << " " << m_Vz[i] << endl;
    file.close();

    //************************************************************************************************************************
    //выводим профили температуры и скорости.
    //************************************************************************************************************************
    file.setFileName(QString("/T(r)_1_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0125));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0125)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/T(r)_2_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0250));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0250)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/T(r)_3_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0375));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0375)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/T(z)_1_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0625,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.0625,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();

    file.setFileName(QString("/T(z)_2_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0750,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.0750,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();

    file.setFileName(QString("/T(z)_3_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0875,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_T, geometry::_2d::Point(0.0875,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();

    file.setFileName(QString("/Vz(r)_1_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0125));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0125)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/Vz(r)_2_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0250));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0250)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/Vz(r)_3_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0375));
        ostream << 0.05 + i*(0.05)/100 << " " <<
                   fem::triangles::functionValue(triangle, m_Net, m_Vz, geometry::_2d::Point(0.05 + i*(0.05)/100,0.0375)) <<
                   endl;
    }
    file.close();

    file.setFileName(QString("/Vr(z)_1_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0625,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(0.0625,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();

    file.setFileName(QString("/Vr(z)_2_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0750,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(0.0750,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();

    file.setFileName(QString("/Vr(z)_3_4.dat").prepend(m_OutputPath));
    file.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < 101; i++)
    {
        int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(0.0875,i*(0.05)/100));
        ostream << fem::triangles::functionValue(triangle, m_Net, m_Vr, geometry::_2d::Point(0.0875,i*(0.05)/100)) <<
                   " " << i*(0.05)/100 << endl;
    }
    file.close();
}
